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Abstract 

The computer-assisted modeling of re-entrant production lines, and, in particu- 
lar, simulation scalability, is attracting a lot of attention due to the importance 
of such lines in semiconductor manufacturing. Re-entrant flows lead to competition 
for processing capacity among the items produced, which significantly impacts their 
throughput time (TPT). Such production models naturally exhibit two time scales: 
a short one, characteristic of single items processed through individual machines, 
and a longer one, characteristic of the response time of the entire factory. Coarse- 
grained partial differential equations for the spatio-temporal evolution of a "phase 
density" were obtained through a kinetic theory approach in Armbruster et al. [2]. 
We take advantage of the time scale separation to directly solve such coarse-grained 
equations, even when we cannot derive them explicitly, through an equation- free 
computational approach. Short bursts of appropriately initialized stochastic fine- 
scale simulation are used to perform coarse projective integration on the phase 
density. The key step in this process is lifting: the construction of fine-scale, discrete 
realizations consistent with a given coarse-grained phase density field. We achieve 
this through computational evaluation of conditional distributions of a "phase ve- 
locity" at the limit of large item influxes. 



Key words: Production line, Re-entrant, Equation-free, Coarse projective 

integration 

PACS: 05.45.-a 



Preprint submitted to Elsevier Science 



2 February 2008 



1 Introduction 



Semiconductor production lines are billion dollar investments and their per- 
formance characteristics are carefully studied and modelled in detail through 
discrete event simulations. In such simulations, a network of processors (ma- 
chines) is set up with connections that describe the product flow through a 
factory. As the production process of a typical semiconductor goes through 
several similar layers, machines may process an item several times at different 
stages in the production process (re-entrancy) . The passage of an item through 
a factory is characterized by its throughput time (TPT), which may vary due 
to several factors, most importantly the stochastic features of the machines 
(their failure statistics) and the interaction with human operators. For a typ- 
ical discrete event model, therefore, the time through a factory process is a 
random variable sampled from a (given) probability distribution. 

A basic model describing the flow of a single type of products through a 
processing factory is given by 

(d) 6 n CL n + T n , 



where a n is the time when item n enters the factory and e n its exit time. 

The time interval r n for which the item stays in the processing factory is sam- 
pled from the distribution T(r,t = a n ), which is determined at the entrance 
time. In engineering practice, the most significant influence on the form of 
this distribution comes from the total number of items in progress, i.e., the 
work in progress (WIP) (see [1,2]). As a result, one rewrites T(r,t = a n ) in 
the form of T(r, WIP), where WIP is the current WIP at the time a n . When 
we fix the TPT at the beginning of the simulation, we essentially treat the 
factory as a single queue whose length at item arrival determines the time 
that the item needs to get processed. However, for re-entrant flows, significant 
changes in WIP during the time interval [a n , e n ] may lead to a change in the 
TPT. In [2], this situation is treated by introducing the concepts of phase (a 
scaled position) and phase velocity. The latter is stochastically updated from 
a distribution that depends on the total WIP in the factory. In this manner, 
the TPT is a variable that can dynamically change, impacted by later arrivals 
of new items. 

In [2] an advection-diffusion equation for the phase density p(x, t) as a function 
of the phase, x, and the time t is derived via a Chapman- Enskog expansion 
of a Boltzmann equation. This PDE can be solved in a straightforward way 
allowing us to determine all relevant quantities: total Work in Progress (WIP), 
the WIP distribution, the throughput and the throughput time. The relevant 
time scale for the PDE simulation is the time scale of the TPT evolution of the 



(b) 



dV{r n < r} = T(r, t = a n )dr. 
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whole factory. In contrast, for a discrete event simulation model, the relevant 
time scale is that of processing by an individual machine, which is at least an 
order of magnitude smaller than the TPT. In addition, in order to obtain a 
meaningful statistical average, a large ensemble of items must be used. As a 
result, the discrete event simulations are time consuming, and the simulation 
effort scales with the number of steps in the factory and the number of items 
produced. In contrast, the continuous model equation is deterministic and its 
simulation is independent of the number of steps and the number of items 
produced. 

Comparison of a discrete event simulation (DES) model detailed below and 
the explicitly derived advection-diffusion equation was performed in [2] and 
the results were in reasonable agreement; in that case, the PDE model was 
derived "from first principles" . However, the particular discrete event simula- 
tions already constituted a significant approximation, as they were written for 
particles in a continuum. In the general case of a large scale DES with many 
items and many machines, global continuum evolution equations for the item 
density cannot easily be analytically derived in a closed form. The question 
then arises whether we can still observe the time evolution at the density 
level, maintaining the scalability advantages of a density-type model. The re- 
cently developed Equation-Free (EF) approach is aimed at precisely this type 
of problem: we have a fine-scale (here DES) simulator, we suspect that we can 
model its coarse-grained behavior at the level of a density evolution equation, 
but do not have this equation explicitly available. The approach attempts to 
solve the equation (integrate it, find its stationary solutions and their sta- 
bility, etc.) by designing short bursts of appropriately initialized simulations 
with the fine-scale (here DES) model. The quantities required for scientific 
computation with the unavailable density-level equation are then estimated 
on demand horn these short fine-scale runs [15,10,9]. In this paper we will use 
a particular equation-free algorithm (coarse projective integration [4,5,14]) to 
illustrate the application of the approach to the re-entrant line problem. 

Other computational tasks like coarse-grained fixed point computation, sta- 
bility, bifurcation analysis, control, and computation of self-similar solutions, 
have been demonstrated in the literature for various types of "inner", fine- 
scale simulators [7,11,12,3,16]. A detailed discussion of the methods can be 
found in [10,9]. A key step in EF computations is the so-called lifting: the 
construction of fine-scale states consistent with given values of their coarse- 
scale observables (the dependent variables in the unavailable coarse equation). 
This is required in order to initiate new short bursts of fine-scale evolution. 
This is not a deterministic step - many fine-scale configurations share the same 
observables. However, the assumption that an evolution equation exists and 
meaningfully closes at the level of these coarse observables suggests that the 
details of various lifting realizations do not affect the long-time coarse-scale 
evolution of the observables (see [7] for detailed discussions). In applying the 
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EF methods to re-entrant factory production, one may use the joint "phase 
and phase velocity" density f(x,v,t) as the observable of choice. However, 
the work in [2] suggests that at an appropriate (large influx) parameter limit, 
one can close an equation in terms of a simpler coarse-grained observable: just 
the phase density p(x,t). While for the former the lifting procedure would 
be straightforward (just generating two random variables according to their 
joint probability density function [16]), for the latter care must be taken in 
generating joint probability densities based on just the phase density p(x,t). 

The paper is organized as follows. In Section 2 we briefly describe the phase 
model of a re-entrant factory discussed in [2]. Section 3 contains our main 
equation-free results: our lifting procedure, as well as coarse projective in- 
tegration based on "finite-difference" -type observations of the phase density 
evolution. We conclude with a brief summary, discussion, and possible exten- 
sions of the approach. 

2 A Phase Model for a Re-entrant Factory 



2. 1 The Discrete Model 



The phase s of an item in a given factory of a supply chain is defined as the 
antiderivative of the phase velocity, l/r(t), where r{t) may change with time, 
- in contrast to the constant r n in Eqn. (1). With a deterministic TPT r(t) 
the improved model reads: 



(a) s = 0(t), 

(6) ^ = i 0K)=O, t>a n . (2) 

The exit time e n is the time at which (p(t) = 1. When r(t) is constant, model 
(2) reduces to model (1). When the TPT is a random variable sampled from a 
distribution T(r,t), a discrete form corresponding to (2) must be used (where 
uj is an update frequency): 



0(t+I) = 0(t) + _ L- 

UJ UJT{t) 

<f)(a n ) = 0, t > a n , dV{r(t) < r} = T (r, t)dr. (3) 

Qualitatively, the value of uj is influenced by the number of items entering the 
factory within a characteristic TPT scale. 

Given a prescribed TPT distribution function T(r, t) (which we assume can 
be expressed in terms of WIP as T(r,WIP)), an executable algorithm for 
evolving the phase s and the real-time TPT r(t) as well as the exit time e n is 
also given in [2], namely, 
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<j>{t + At) = 0(t) 



r(ty 

r(t + At) = K(t)r)(t) + (1 - «(t))r(t), t > a n , 

(b) V{K(t) = 1} = ouAt, V{K{t) = 0} = 1 - cuAt, 
dV{rj(t) < r} =T(r,t)dr, 

(c) 0(a n ) = 0, dV{r(a n ) <r} = T (r, a n )dr, (4) 

where the update frequency ui depends on the throughput time r(t) and the 
time t itself, i.e., u = uj(r(t),t). The simulation time step At is chosen such 
that At < -. For small enough At this algorithm can be used to numerically 
solve model (3). 

In [2], u is chosen as 

-M) = ^ (5) 



where \(t) is the influx of items and T_i = Jr 1 T(r,t)dr. 

The procedure to compute the phase and throughput time (TPT) of an item 
in the factory is then as follows: 

(1) Set the item's initial phase to as it enters the factory; increase the 
WIP W(t) by one and adjust the distribution of TPT, T(r, t) accordingly 
(remember that we have assumed that this distribution can be written 
as T(r, WIP)); the item's initial TPT is then sampled from the udpated 
T(r, WIP); 

(2) Compute cj(r(t),t) as above, and sample the parameter n(t) according 
to its distribution (4)(b); 

(3) Compute the TPT r{t + At) and the phase (fi(t + At) according to (4) (a); 

(4) Update the distribution of TPT, T(r, t), according to current WIP at the 
time t + At, and go back to Step (2); repeat this loop until 0=1. The 
time when = 1 is the exit time of the item, e n . 

2.2 Density Equations 

Let the joint number density of phase and phase velocity, f(x, r, t), be defined 
as f(x,r,t) = 9 F ^(tyQr~ r '^ ; where F(<f> < x,t < r,t) is the number of items 
whose phase < x and TPT r < r. The PDE for f(x, r, t) is derived in [2]. It 
is given by 



df ldf 

+ = T(r, t) I u(r', t)f(x, r', t)dr' - w{r, t)f(x, r, t), 
x>0, t > 
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f(0,r,t)=rX(t)T(r,t), /(z,r,0) = 



(6) 



weakly in x, r and t, where X(t) is the influx of items. The chosen form of 
uj in (5) guarantees that the discrete model (1), when there are no temporal 
changes in influx and TPT, is a Monte Carlo scheme corresponding to a par- 
ticular solution of (6), and that the instances of random number generation 
are roughly the same for models (1) and (4). 

If we define the number density of phase, p(x,t), as p(x,t) = J f(x,r,t)dr, 
then in the limit that XqCTj- ^> 1 (where Ao and a\ are characteristic scales 
of influx and standard deviation of TPT, respectively) one can write a closed 
equation for its evolution as follows: 



dp dF 
dt dx 



0, x, t > 0, 

1 T_! 1 d(T 2 /Ti] 
T x X T x dt 
F(0,t) = X(t), p{x,0) = 



F(x,t) = C(t)p-D(t) 



dp 

dx' 



C 



D 



T-i T 2 - T[ 
X T? 



(7) 



The moments of the probability distribution T(r, t) appearing in this formula, 
Ti(t),i= -1,1,2, are defined by T^t) = J r l T (r, t)dr, i = -1,1,2. 

Based on the proof of Theorem T3 in [2], it can also be deduced that in this 
limit, the joint number density f(x,r,t) is controlled by (slaved to) just the 
number density of phase, p(x,t). Their relationship is given by 

K, T , t) = e2g£ ( 8) 

Equation (8) gives a simple relationship between the joint number density and 
the number density of phase. It implies that the conditional number density 
f(r\x,t) (= f(x,r,t)/p(x,t)) is a constant rT(r,t)/Ti with respect to the 
phase coordinate x. 

3 Equation-Free Analysis for the Two-scale System 

The Equation-Free approach consists of an ensemble of computational tools 
that can be used to study the coarse-grained, macroscopic behavior of sys- 
tems based on their underlying fine-scale, microscopic simulators (see e.g. 
[4,7,11,11,3,16]). The basic element of these algorithms is the so-called coarse 
time-stepper, whose role is to connect observables across different scales. The 
coarse time-stepper consists essentially of three components: lifting and micro- 
simulation followed by restriction. The lifting is a procedure that generates 
micro-scale realizations of a system state consistent with given values of their 
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macro-scale observables; while restriction is the reverse: obtaining macro-scale 
observables from the fine-scale system state. Since a fine-scale system state 
normally possesses far more degrees of freedom than a few macro-scale ob- 
servables, the lifting procedure is not a one-to-one mapping in general. Care 
needs to be taken when a lifting algorithm is implemented and tests are usu- 
ally required to check if the macro-scale evolution is sensitive to the details of 
a particular lifting. 

If we have reason to believe that useful macroscopic equations accurately close 
at the level of a few macroscopic observables (e.g. in terms of a few moments of 
microscopically evolving distributions, as is the case in hydrodynamic equation 
derivation from Boltzmann- level descriptions), equation-free schemes allow us 
to solve the coarse-grained equations without the explicit closures needed to 
write them in closed form. The idea is to use short bursts of fine-scale simu- 
lation to evaluate the right hand sides of the unavailable closed equations on 
demand. Traditional continuum numerical analysis is thus transformed into 
protocols for the design and processing of repeated short bursts of compu- 
tational experiments with the fine-scale solver. The assumption that coarse- 
grained equations in principle exist and close in terms of a few coarse-grained 
observables implicitly suggests that the details of a lifting should be quickly 
forgotten; the role of the brief fine-scale simulation is to "implement" this loss 
of memory of initial features of the fine-scale state. 

3.1 The Lifting Step in the Equation-Free Approach 

We are interested in enabling short bursts of simulation with the discrete 
phase model (4) to numerically analyze the evolution of the number density 
p(x,t) - a coarser-grained observable than the full f(x,r,t). This should only 
be attempted at conditions when a deterministic equation closes with p(x,t), 
i.e. when the evolution of the full f(x,r,t) is controlled by p(x,t). For such 
a closed equation to exist, it should be possible (possibly after a short initial 
transient) to express f(x,r,t) in terms of p(x,t). This section investigates the 
effect of the parameter Ao<Xj- on the relationship between p(x, t) and f(x, r, t). 
This is both because A <Xj- is the only dimensionless parameter that appears 
in the discrete model, and also because we are interested in cases of large item 
influx within the characteristic TPT scale, when \qo~^- 3> 1. 

In what follows, the discrete model (4) is executed with constant influxes and 
uniform TPT distributions given in Table 1. For all cases, an ensemble of 
7, 000 realizations is used in order to obtain smooth evolution of the densities 
f(x, r, t) and p(x, t). At time t=16sec, the phase <fi and TPT r of items existing 
in the factory (i.e., < <fi < 1) are recorded and used to plot a surface for the 
conditional density f(r\x,t) (Fig. 1). It is found that when \qo~j- is sufficiently 
large (as in Case 3,5,6,8,9), f(r\x,t) has become independent of the phase 
coordinate x. Actually, for our particular choice of T(r, t) shown in cases 
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Case Number influx, A(t) PDF of TPT, T(r, t) 



1 

1 


U.O 


uniform in [0.1,2] 


n 
2 


10 


uniform in [0.1,2] 


Q 



on 


uniform in [0.1,2] 


1 

4 


U.O 


uniform in [0.1,4] 





1 n 
1U 


uniform in [0.1,4] 


6 


20 


uniform in [0.1,4] 


7 


0.5 


uniform in [0.1,8] 


8 


10 


uniform in [0.1,8] 


9 


20 


uniform in [0.1,8] 



Table 1 Cases of influxes and uniform distributions of TPT both constant over the 
time domain 

3,5,6,8,9 of Table 1, we find that the conditional density f(r\x,t) follows a 
linear distribution in r passing through the origin (Fig. 2). Clearly, in order to 
lift effectively, one needs to know the support of the conditional distribution; 
for our choice of T(r, t) the lower limit of its support in r does not depend on 
WIP or time, and it makes sense to take this also as the lower limit of support 
of f(r\x,t); we have observed that - to within acceptable error - the upper 
limit of support of f(r\x,t) coincided in our simulations with the upper limit 
of support of T(r,t). In this case, f(r\x,t) can be approximated in the form 
of rT(r,t)/C where C is a constant independent of r. Obviously, C equals 
Ti(t) since / f(r\x,t) = 1. Other choices of X(t) and T(r,t) (e.g., a linearly 
increasing influx and a linear TPT distribution) give rise to more or less the 
same relationship between p(x,t) and f(x,r,t). The "constitutive equation" 
(8) and the parameter regime under which it is justified were thus found using 
only short simulations with the discrete fine-level model. 

For a prescribed T(r,t), and starting at a given p(x,t), the lifting algorithm 
can now be formulated when Xqct^- ^> 1 as follows: 

(1) Calculate the WIP(t): WIP = $ p(x,t)dx. 

(2) Calculate the probability density function (or normalized number den- 
sity) f p (x,t) of the phase coordinate: f p (x,t) = ^jp- 

(3) Since the number of items generated can only be an integer, we have to 
systematically generate an ensemble of integers whose mean value equals 
WIP. Let a = int(WIP) + 1 - WIP, where int(WIP) is the maxi- 
mum integer not greater than WIP. Select a random variable p which 
is uniformly distributed in [0, 1]. If p < a, then the number of items is 
int(WIP). Otherwise it is int{WIP) + 1. 

(4) Compute the cumulative distribution function (CDF) of the phase coor- 
dinate Fx(x,t) = Jo f p (x,t) and its inverse IFx(Fx,t). 
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Case 2 



Case 3 












Fig. 1. Conditional distributions of TPT, f(r\x,t) 
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Fig. 2. Conditional TPT distribution, f(r\x,t = 16s). Top: Case 5; Bottom: Case 9 



(5) For each realization in the integer ensemble, generate the phase coor- 
dinates of items: (pi = IFx(0,t),i = 1,2, • • ■ ,n, where n is the integer 
number of items and d, i = 1, 2, • • • , n are random numbers uniformly 
distributed in [0,1]. 

(6) Compute the CDF of the TPT F R (r,t) = JI^ f(r\x,t) and its inverse 
IF R (F R ,t). The phase coordinate x is suppressed in F R (r,t) since it is 
independent of x. 

(7) For items in each realization, generate their respective TPT's: n = IF R (ipi, t), 
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i = 1, 2, • • • , n, where ipi are also random numbers uniformly distributed 
in [0,1]. 

The following test was performed to validate the above procedure. A true 
trajectory of p(x, t) for Case 9 in Table 1 is computed using 5, 000 ensemble 
realizations in the time domain directly by the discrete phase model (4). Then 
the trajectory is interrupted at t — 10; the number density p(x,t) at the 
interruption is retained, and then lifted as described above to start a new 
discrete evolution. We observed that the relative difference between the two 
trajectories (the normally continued one and the one starting from the lifting 
after the interruption) is quite small (within 2%) even just after the lifting. 
Another realization of the lifting may be taken which uses a different random 
number seed. The result shows no significant difference than the previous 
situation. It can therefore be concluded that the lifting algorithm is effective 
and particular liftings do not significantly affect subsequent restricted number 
densities. 

3.2 Coarse Projective Integration of the Density-level System 

Coarse Projective Integration (CPI) is a numerical technique developed in the 
EF framework to evolve in time the coarse-grained observables of a multiscale 
system, estimating their temporal derivatives using the underlying fine-scale 
simulator [4,5,14]. This technique is suitable for systems whose coarse- and 
fine-level temporal scales are well separated, i.e., the coarse-level observables 
are smooth over a temporal scale that is significantly larger than the fine-level 
evolution scales (the scales that it takes for higher-order system observables 
to become slaved to the slow, "master" ones). In the same spirit with adaptive 
time-step selection based on error control for standard deterministic integra- 
tors [13], real-time computational tests for projective step selection have to be 
taken to ensure error control in the CPI evolution of macroscopic observables. 

We now implement CPI for our two-scale supply chain model. The prescribed 
influx and TPT distribution are shown in Fig. 3, respectively. The time step 
Atj in the discrete model is chosen as Std = 10~ 3 sec. In the following, a finite- 
difference motivated representation is used to evolve coarse-grained observ- 
ables, p(x,t), through CPI. If an explicit equation for p(x,t) was discretized 
in space with finite differences, we would evolve, in time, the values of the field 
at a number of points - and we would use differences between these values 
at each moment in time to approximate the right-hand side of the explicit 
equation. Here, we do not use the finite point representation to approximately 
evaluate the right-hand side of the equation - we use the representation in- 
stead to lift to an item distribution in the entire domain, and run the fine-scale 
simulator to estimate the left-hand side of the evolution of the representation. 
It is these values that integration codes need to solve the initial value prob- 
lem, and the way they process the numbers does not depend on whether they 
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t WIP 



Fig. 3. Influx and probability distribution of the throughtput time 

came from approximating the right-hand side of the equation, or from esti- 
mating its left-hand side. We projectively integrate the number density p(x, t) 
at equally spaced phase points Xj = j/M,j = 0,1,- • • , M; here we choose 
M = 8. (We have also successfully performed coarse projective integration us- 
ing as variables the projections of the solution on the first few of its empirical 
orthogonal global basis functions (POD modes) in a POD-assisted approach; 
these results are not shown because of space limitations). At the starting time 
t , the number density p(xj, t ),j = 0, 1, • • • , 8 is lifted (according to the lifting 
algorithm in Section 3.1) to initiate the discrete phase model (4). This model 
is subsequently evolved for 20 discrete time steps. At steps i, i — 12, 14, • • ■ , 20, 
phases of items in progress are restricted and binned to obtain number density 
histories p(x,t + iStd). These number density histories are used to approxi- 
mate the temporal derivatives dp(x,t)/dt at time t + 205td via least-squares 
fitting. The number density after a coarse-level time interval At c can then be 
obtained based on a simple forward Euler explicit algorithm by projection: 

p(x, t + At c ) = p(x, to + 205t d ) + (At c - 205 d )^- {x 



dt 

This explicit Euler scheme illustrates the simplest way to project the number 
density over a coarse time interval, here with computational savings of roughly 
(1 — 205d/A c ). More sophisticated projective algorithms, including projection 
templated on higher-order continuum integration schemes [14], taking advan- 
tage of multiple spectral gaps [6] and even coarse implicit schemes [7], can be 
utilized as well. 

Figure 4 compares trajectories of WIP and outflux in true and CPI evolutions 
for some selections of the coarse time step At c . All computations involve an 
ensemble of 5, 000 realizations. It is found that the CPI trajectories of WIP 
almost coincide as At c = 0.2s and At c = 0.3s, which suggests that the result 
at At c = 0.2s is correct. Comparison also displays that WIP and outflux at 
At c = 0.2s indeed match their true trajectories. Number densities of the true 
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time time 

Fig. 4. True trajectories and some CPI evolutions of WIP and ouflux. Left: WIP, 
clusters of points represent WIP's at time steps ikAtj, i = 0, 1, • ■ ■ , 10 immediately 
following the lifting; Right: outflux. 
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Fig. 5. Number densities of the true evolution and CPI evolution with At c = 0.2s. 
Left: true evolution; Right: CPI result 

evolution and the CPI result with At c = 0.2s are also shown in Fig. 5. By using 
coarse projective integration, the computational load is substantially reduced 
by 78.3%, taking into account the trial efforts in selecting an appropriate 
coarse timestep size for projection. 

4 Conclusions 

We demonstrated certain features of Equation-Free coarse-grained computa- 
tion for a re-entrant supply-chain model. In cases where an explicit equation 
for the processed item number density cannot be easily derived, simulations 
with the fine-scale model can be used to determine the conditions under which 
such an equation can, in principle, exist. Once these conditions are established, 
equation-free methods such as coarse projective integration can be used to 
evolve the coarse-level density through short bursts of appropriately initial- 
ized runs with the fine-scale simulator. The lifting process essential to this 
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EF approach was investigated and a particular governing form establishing a 
relation between the joint number density f(x,r,t) and the number density 
p(x, t) was found; the form is of course only valid for this model. 

EF algorithms based on matrix-free iterative linear algebra (like Newton- 
Krylov GMRES, [8]) can also be used to effectively implement contraction 
mappings to find stationary item number densities, if the influx to a factory 
reaches a stationary state at long times. Additional tasks like continuation, sta- 
bility and parametric sensitivity analysis, and even control and optimization 
computations can in principle be implemented in an equation-free framework. 

Our illustrative example involved a Monte-Carlo type solution of a Boltzmann 
equation, which, at a certain limit, approaches a discrete event simulator. In 
this case, both the Boltzmann equation and the reduced equation for the item 
number density at the appropriate limit were analytically available; this al- 
lowed us to validate our equation-free computations. The real challenge lies in 
wrapping equation-free algorithms around true discrete event simulations for 
realistic processing factory configurations. While the derivation of continuum- 
level equations may be extremely difficult or practically infeasible in such 
cases, our computer-assisted approach remains essentially the same, indepen- 
dent of the details in the underlying fine-scale simulator. We believe that this 
approach holds promise for facilitating the extraction of system-level informa- 
tion from complex discrete event simulators. 
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